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1  INTRODUCTION 

This  document  introduces  a  user  interface  that  provides  fast  and  efficient  general  least- 
squares  analysis  of  experimental  data.  It  is  suitable  for  use  in  the  Sunwindows  environment  on 
a  Sun  UNIX  workstation.  The  package  is  capable  of  providing  solutions  to  both  linear  and 
non-linear  least-squares  models  of  experimental  data,  and  incorporates  full  graphics  and 
error  analysis.  Known  uncertainties  in  the  data  are  easily  incorporated  as  "weights”  for  each 
experimental  point  if  required. 

The  routine  is  user  specific  and  as  such  a  limited  knowledge  of  UNIX  as  well  as  C  and/or 
Pascal  is  necessary.  However  an  attempt  has  been  made  to  lessen  this  burden  as  much  as 
possible  without  destroying  the  flexibility  of  the  program.  An  overview  of  the  operational 
details  of  the  package  is  given  in  Section  2.  Examples  of  how  to  use  the  package  are  given  in 
Section  3. 

The  linear  least-squares  analysis  allows  unweighted,  y-weighted  and  x-  and  y-weighted 
optimisation,  and  in  the  former  cases  either  the  slope  or  intercept  may  be  fixed  if  required. 
Since  the  theory  of  unweighted  and  y-weighted  linear  least-squares  is  well  known[l],[2],  it 
will  not  be  described  in  this  document.  The  theory  for  the  case  where  both  the  dependent  and 
independent  variables  exhibit  uncertainties  is  not  well  known  however,  and  as  such  is 
outlined  in  Appendix  A. 

Appendix  B  details  the  theory  of  non-linear  least-squares  optimisation  using  the  Levenberg- 
Marquardt  algorithm,  which  is  the  routine  implemented  in  this  package.  The  routine 
incorporates  options  for  unweighted  or  y-weighted  fitting  of  models  with  multiple 
parameters  (maximum  specified  by  the  user),  any  of  which  may  be  easily  fixed  or  floated  as 
desired.  An  outline  of  the  method  used  to  estimate  the  errors  in  the  final  parameters  is  also 
given. 

The  package  is  public  domain  software  and  as  such  there  is  no  guarantee  it  is  bug  free.  It  may 
be  freely  distributed  and  augmented,  the  only  request  being  that  the  original  and  subsequent 
substantial  authors'  names  remain  at  the  top  of  the  major  routines. 
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2  OPERATIONAL  DETAILS 

In  this  section  a  description  of  the  system  and  the  details  of  establishing  an  appropriate 
environment  are  given.  Figure  1  shows  an  overview  of  the  environment  in  which  spg_fit 
operates.  The  role  of  each  of  the  modules  is  discussed  below. 


USER 

INPUT 


DATA  PARAMETER  PLOT 

FILES  FILES  FILES 


Figure  1  The  file  environment  of  spg_fit 
2.1  The  function  file 

The  function  file  supplies  both  the  model  and,  optionally,  the  derivatives  of  the  model 
with  respect  to  each  parameter.  The  function  may  be  written  in  either  C  or  Pascal. 
Note  however,  that  if  the  routine  is  written  in  Pascal,  all  call  parameter  arrays  must 
be  var  since  the  main  spg_fit  C  routine  communicates  using  pointers.  Details  of  the  call 
parameters  may  be  found  in  README  included  on  the  archive  and  given  in  Appendix  C. 
The  filenames  must  be  either  cfn.c  or  pfn.p  for  the  C  and  Pascal  routines  respectively. 
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Examples  of  both  C  and  Pascal  routines  are  included  on  the  tape  (in  the  directory 
Fit_lib)  and  are  given  in  Appendix  D.  They  include  simple  models  such  as 
trigonometric,  exponential  and  distribution  functions.  The  user  must  copy  the  required 
routine  into  the  appropriate  file  (cfn.c  or  pfn.p)  and  recompile  (see  Section  2.7).  If  the 
spg_fit  routine  is  recompiled  without  the  presence  of  the  function  file  in  the  current 
working  directory,  a  dummy  function  will  be  automatically  linked.  In  this  case  only 
the  linear  fitting  mode  will  be  operational  (see  Section  2.6.2). 

As  noted,  supplying  the  first  derivatives  of  the  function  with  respect  to  each  parameter 
is  optional.  If  the  analytic  derivatives  are  supplied,  the  program  will  run  slightly 
faster  and  the  usual  round-off  errors  associated  with  numerical  evaluation  of  the 
derivatives  will  be  avoided.  However,  if  the  derivatives  are  not  supplied  in  the 
function  file,  spg_fit  will  automatically  recognise  this  and  numerical  evaluation  of  the 
derivatives  will  be  undertaken. 

2.2  The  data  file 

The  data  file  is  a  text  file  containing  the  experimental  data  and  optionally  any 
experimental  variances.  The  format  of  the  data  file  is  simply: 

2  2 


XI 

yi 

CT*1 

CTyi 

2 

2 

x2 

y2 

ax2 

CTy2 

xn 


yn 


If  there  are  only  three  columns,  the  third  column  is  automatically  interpreted  as  a 


and  if  there  are  only  two  columns  the  program  recognises  that  the  uncertainties  of  the 
individual  data  points  are  unknown.  The  columns  may  be  separated  by  either  spaces  or 
tabs.  For  graphics  purposes,  the  data  (x)  should  be  presorted  into  either  ascending  or 
descending  order. 


2.3  The  function  information  file 

The  function  information  file,  called  fn_info,  is  a  text  file  that  contains  the  names  of 
the  model  parameters,  entered  line  by  line.  The  routine  uses  the  text  strings  from  this 
file  for  display  purposes  and  also  to  calculate  the  number  of  parameters  in  the  model. 
It  is  only  required  for  the  non-linear  mode  of  operation.  It  is  important  that  the  order  of 
the  parameter  names  be  the  same  as  the  order  in  which  the  parameters  are  specified  in 
the  user's  routine.  Similarly,  it  is  important  to  remember  to  change  the  fn_info  file 
whenever  the  non-linear  function  is  changed. 
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2.4  The  curve  save  file 

The  curve  save  file  is  named  and  established  interactively  from  the  main  control  panel 
and  contains  the  theoretical  curve  whenever  "Save"  is  pressed.  Saving  the  curve  is  not 
done  automatically. 

2.5  The  parameter  save  file 

The  parameter  save  file  is  again  named  and  established  interactively  from  the  main 
control  panel.  This  feature  is  useful  not  only  for  storing  the  final  analysis  results  but  can 
be  used  to  provide  an  initial  guess  for  parameter  values  based  on  a  previous  fit.  This  is 
achieved  by  specifying  the  appropriate  parameter  save  file  on  the  main  control  panel 
and  clicking  the  "Load"  button.  When  "Save"  is  clicked,  the  specified  file  will  contain 
the  parameters  as  well  as  their  corresponding  errors. 

2.6  The  interface 

The  interface  consists  of  a  main  control  window  and  display  canvas,  a  "linear  fit"  panel 
and  a  "non-linear  fit"  panel.  The  "fit"  panels  are  pop-up  windows  selected  by  choosing 
either  the  "linear"  or  "non-linear"  options  on  the  main  control  panel  via  the  "select 
model  type"  toggle.  The  function  and  layout  of  each  of  the  three  main  panels  is 
described  in  the  following  sections. 

2.6.1  The  main  control  window  and  display  canvas 

The  main  control  window  and  display  canvas  appear  immediately  the 
programme  is  executed.  Figure  2  shows  the  layout  of  this  window  along  with 
the  linear  fit  options  panel. 

The  field  Directory  contains  the  current  working  directory  which  may  be 
changed  at  any  time.  Files  for  data  and  saving  may  be  entered  into  the 
indicated  fields.  The  toggle  shown  in  the  top  centre  allows  switching 
between  the  non-linear  and  the  linear  modes.  Figure  3  shows  the  non-linear 
selection. 

The  display  graphics  has  two  modes  of  operation: 

•  If  a  datafile  is  specified,  the  routine  will  calculate  ranges  and 
display  the  data  points  with  any  associated  error  bars  (the  error  bar 
length  is  defaulted  to  two  standard  deviations  either  side  of  the 
point).  Once  a  parameter  has  been  specified  a  curve  will  be  drawn 
spanning  the  data. 

•  Alternatively  if  no  data  is  specified  but  ranges  and  parameters  are 
correctly  input,  the  routine  will  display  the  model  curve  over  the 
selected  horizontal  range  of  the  graph  on  the  display  canvas.  This  is 
useful  when  examining  the  behaviour  of  the  model. 


4 


UNCLASSIFIED 


UNCLASSIFIED 


ERL-0544-GD 


UNCLASSIFIED 


6 


UNCLASSIFIED 


UNCLASSIFIED 


ERL-0544-GD 


The  Curve  Points  prompt  on  the  side  panel  refers  to  the  number  of  points  the 
routine  will  use  to  draw  the  curve.  The  default  is  set  to  100,  or,  when  there  is 
data,  to  the  number  of  input  data  points.  Problems  with  aliasing  will  only 
affect  the  display  and  not  the  fitting.  During  the  optimisation  this  number 
will  be  equal  to  the  number  of  data  points,  since  spg_fit  will  obtain  the  values 
of  these  points  from  the  Levenberg-Marquardt  call.  Upon  completion  the 
number  of  curve  points  will  revert  to  that  specified. 

2.6.2  The  linear  fit  panel 

The  linear  fit  panel  will  appear  upon  selection  of  the  "linear  fit"  option  from 
the  toggle  on  the  top  control  panel. 

Figure  2  shows  the  linear  least-squares  mode  of  operation.  Beside  each 
parameter  is  a  box  which  may  be  checked  or  left  clear  if  floating  or  fixing  of 
the  parameters  is  required.  A  weighted  correlation  coefficient  (R)  is 
calculated  and  displayed,  as  well  as  a  value  for  the  (reduced)  weighted  sum 
of  the  square  of  the  residuals  (reduced  chi  squared).  The  toggle  has  three 
choices,  unweighted,  y-weighted  and  x-  and  y-weighted.  Only  in  the  latter 
case  is  the  option  to  fix  parameters  not  included.  By  fixing  the  slope  at  zero, 
this  panel  may  be  used  to  calculate  either  the  weighted  or  unweighted  mean 
of  the  data,  along  with  the  associated  uncertainty. 

2.6.3  The  non-linear  fit  panel 

The  non-linear  fit  panel  will  appear  upon  selection  of  the  "non-linear  fit" 
option  from  the  toggle  on  the  main  control  panel.  This  is  shown  in  Figure  3. 
This  mode  of  operation  may  be  used  when  the  model  is  non-linear  with 
respect  to  the  parameters  or  for  generalised  linear  least-squares  analysis. 

The  layout  and  operation  is  similar  to  the  linear  case.  The  values  of  the 
parameters  may  be  changed  at  any  time,  with  the  routine  updating  the  screen 
automatically.  The  y}  value  corresponding  to  the  new  parameter  value  will 
also  be  displayed  if  the  Curve  Points  (as  described  in  Section  2.6.1)  is  equal  to 
the  number  of  data  points.  This  is  so  that  two  function  calls  are  not  required, 
ie.,  one  for  the  plot  and  one  for  the  yp-  calculation,  as  this  would  be  time 
consuming  for  complex  functions.  If  Curve  Points  is  not  equal  to  the  number  of 
data  points,  then  the  y}  value  is  only  updated  once  the  non-linear  fitting  is  in 
progress. 

Once  the  floating  parameters  are  chosen  and  the  weighting  option  selected, 
"Fit  Selected"  may  be  chosen  to  begin  the  optimisation.  "Fit  All"  will  select 
all  parameters  to  float.  If,  after  an  iteration,  the  fit  is  found  to  have 
improved,  the  routine  will  update  the  parameters  on  the  screen,  including  the 
Levenberg  Marquardt  parameter  (lambda)  and  the  value  for  y}.  By  pressing 
the  LI  (Stop)  key  on  the  Sun  keyboard  the  fit  may  be  halted  at  any  stage. 
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This  is  not  immediate  since  the  routine  waits  for  any  internal  user  function 
calls  to  be  completed.  The  parameter  lambda  may  be  specified  before  fitting 
to  optimise  the  convergence,  however  this  should  rarely  be  necessary. 

Once  optimisation  is  completed  the  "Error  Analysis”  button  may  be  chosen  to 
calculate  the  uncertainty  in  each  parameter.  This  is  taken  as  two  standard 
deviations  (95%),  however  this  may  be  altered  by  changing  the  appropriate 
field  in  "spg_fit.h". 

2.7  Setting  up  the  spg_fit  environment 

The  user  must  execute  a  small  number  of  steps  in  order  to  customise  the  routine  for 
his/her  purposes.  They  are: 

1.  Create  the  source  file  (cfn.c  or  pfn.p),  discussed  in  Section  2.1,  containing  the 
non-linear  model.  This  is  not  necessary  when  a  linear  least  squares  fit  is  all 
that  is  required. 

2.  Create  the  corresponding  fn_info  discussed  in  Section  2.3.  If  this  does  not  exist 
the  routine  will  continue,  however  only  the  linear  least  squares  mode  will  be 
operational. 

3.  Customise  the  compile  parameters.  These  parameters  are  found  in  the  files 
Makefile  and  spg_fit.h  (see  Appendix  C)  and  should  be  checked  before 
compiling. 

4  Re-compile  using  a  simple  make  command. 

Details  of  all  these  steps  may  be  found  in  the  README  file  supplied  on  the  archive, 
and  in  Appendix  C. 

2.8  The  algorithms 

It  is  the  intention  of  this  package  to  provide  a  convenient  means  of  undertaking  the  task 
of  general  least-squares  modelling  in  as  simple  a  manner  as  possible  for  the  user.  It  is 
also  the  intention  of  the  authors  that  the  code  be  accessible  to  the  user  if  changes  or 
updates  are  considered  necessary.  It  is  useful  therefore  to  have  some  understanding  as 
to  how  the  optimisation  routines  used  in  spg_fit  actually  work.  To  this  end,  in  order  to 
make  this  document  as  self-contained  as  possible,  Appendices  A  and  B  contain  a  basic 
outline  of: 

(i)  linear  least-squares  optimisation  for  the  case  when  both  the  dependent  and 
the  independent  variables  are  subject  to  some  uncertainty. 
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(ii)  the  theoretical  basis  of  the  Levenberg-Marquardt  non-linear  optimisation 
algorithm  used  in  spg_fit,  and 

(iii)  a  discussion  of  the  estimation  of  the  statistical  errors  associated  with  the 
optimised  parameters  for  the  case  of  non-linear  functions. 

As  noted  in  the  introduction,  the  common  cases  of  unweighted  linear  least-squares 
fitting  and  of  y-weighted  linear  least-squares  fitting  are  well  known,  and  as  such  will 
not  be  described  in  this  report.  Those  unfamiliar  with  the  basis  of  these  algorithms 
are  referred  to  references  [1]  and  [2],  or  to  any  of  the  standard  undergraduate 
mathematical  texts  dealing  with  numerical  techniques. 

3  AN  EXAMPLE  -  FITTING  A  POWER  FUNCTION  (y=A+B*x") 

In  this  section  an  example  is  followed  through  to  demonstrate  the  procedures  of  using  spg_fit. 

The  first  step  is  to  copy  the  user  function  and  information  files  into  the  appropriate 
filenames.  For  this  example: 

cp  /source_path/Fit_lib/fit_power_fn.c  cfn.c 

cp  /source_path/Fit_lib/info_power_fn  fn_info 

where  source_path  is  the  spg_fit  source  directory. 

Now  the  routine  is  ready  for  compiling.  This  may  be  done  from  any  directory  which  includes 
the  above  two  files.  To  compile  from  a  directory  other  than  the  spg_fit  source  type 

make  -f  /source_path/Makefile  fitc 

If  you  are  working  from  the  source  directory  then  simply  type 

make  fitc 

(fitc  is  used  because  in  this  example  the  function  is  written  in  C). 

This  process  will  leave  the  spg_fit  object  files  in  the  source  directory  and  place  the  function 
object  file  in  the  current  directory.  The  executable  will  reside  in  the  path  specified  by  DEST 
which  is  defined  in  the  Makefile  or  overridden  on  the  command  line. 

The  routine  is  now  executed  by  typing  fit  at  which  time  the  main  control  window  appears.  By 
pressing  over  the  top  centre  toggle  (labelled  "Select  Model  Type” )  with  the  left  mouse  button, 
a  non-linear  options  panel  will  appear  as  shown  in  Figure  4.  From  this  point  the  user  may 
examine  the  properties  of  the  function  or  input  a  datafile  in  the  appropriate  field. 
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In  Figure  5  the  file  "power_noisy"  has  been  entered  and  displayed.  (No  error  bars  for  this 
data).  Note  that  the  number  of  samples  is  automatically  set  to  the  number  of  data  points.  At 
this  point  parameter  values  may  be  changed  manually  to  experiment  with  the  data,  as  shown 
in  Figure  6. 


Figure  5  Sample  Data 
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It  is  initially  desired  to  fit  "power_noisy”  with  all  the  parameters,  hence  either  all  boxes 
are  ticked  and  "Fit  Selected”  chosen  ,  or  simply  the  button  "Fit  All"  is  selected.  With  some 
functions  the  routine  may  fail  if  some  of  the  parameters  are  set  to  zero  at  the  start.  In  the 
power  function  for  instance,  setting  "scale"  to  zero  and  fitting  will  return  IEEE  NaN's.  By 
starting  the  optimisation  with  the  parameter  scale  =  0.1  and  performing  an  error  analysis 
produces  the  screen  shown  in  Figure  7.  The  data  in  this  example  was  generated  with  the 
function  y=x156  with  additive  gaussian  noise  of  variance  0.1.  Note  the  small  value  of  lambda 
at  the  completion  of  the  fit.  Figure  8  demonstrates  a  fit  when  one  of  the  parameters  is 
required  to  be  fixed,  in  this  case  the  intercept.  A  zero  error  is  returned  for  the  fixed 
parameter.  Once  the  fit  is  satisfactory  the  number  of  samples  may  be  increased  and  the  curve 
saved  so  that  a  quality  graphics  package  may  be  used  to  obtain  a  printout. 

The  power  function  example  is  particularly  interesting  in  that  it  also  demonstrates  a  pitfall 
that  is  occasionally  encountered  when  using  the  Levenberg-Marquardt  algorithm.  If  the 
optimised  parameters  shown  in  Figure  7  are  changed  very  slightly  -  say  in  the  last  one  or  two 
digits  -  and  the  fit  re-started,  very  often  the  Marquardt  parameter  will  be  seen  to  "explode", 
there  being  however  little  or  no  associated  change  in  the  parameter  values.  This  is  due  to  the 
fact  that  the  routine  can  have  a  tendency  to  wander  around  the  minimum  if  the  x^ 
hypersurface  exhibits  a  generally  flat  valley  with  complicated  topology  in  the  region  of  the 
optimal  parameter  values[2].  In  this  case  a  very  small  change  in  the  parameter  values  can 
produce  a  large  correction,  which  will  be  rejected  and  the  value  of  X  subsequently  increased. 
This  procedure  can  occasionally  "run  away".  This  is  one  example  where  the  parameter  values 
are  valid  despite  an  associated  large  value  for  the  Marquardt  parameter. 

Due  to  the  possibility  of  the  Marquardt  parameter  "running  away",  a  flexible  ceiling  has 
been  placed  on  the  value  of  the  maximum  number  of  iterations  (see  MAX_ITERATIONS  in 
spg_fit.h).  The  problem  may  be  overcome  by  choosing  different  starting  values,  or  changing 
the  starting  value  of  X. 

Another  problem  the  user  may  encounter  is  the  convergence  of  the  parameters  to  a  local 
minimum.  This  routine  does  not  supply  any  algorithms  such  as  Monte  Carlo  searching  to 
address  this  problem.  Hence  there  must  be  some  judgement  in  choosing  appropriate  starting 
values  for  the  parameters  if  there  are  multiple  minima  in  the  hypersurface.  If  there  are  no 
a  priori  reasons  for  choosing  starting  values  for  a  particular  parameter  then  a  manual  search 
for  other  minima  could  be  undertaken,  taking  note  of  the  values  in  turn. 
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LINEAR  LEAST-SQUARES  FITTING  TO  DATA  EXHIBITING 
UNCERTAINTIES  IN  BOTH  THE  DEPENDENT  AND  INDEPENDENT 

VARIABLES 


A.l  Introduction 

The  solution  to  the  problem  of  finding  the  best  linear  least-squares  fit  to  data  which 
exhibits  uncertainty  in  the  dependent  variable  is  well  known,  and  will  not  be 
considered  here.  The  more  interesting,  and  probably  more  realistic  case  of  finding  the 
"best"  description  of  data  when  there  exists  uncertainty  in  both  the  dependent  and 
independent  variables  is  the  subject  of  this  section. 

The  case  where  both  variables  exhibit  some  uncertainty  has  generally  led  to  the  use  of 
some  approximation  method  to  determine  the  best  fit  to  the  data.  For  example,  one 
may  associate  all  of  the  error  with  the  y  co-ordinate,  calculate  the  usual  weighted 
least-squares  solution  and  then  do  a  similar  calculation  with  all  the  error  associated 
with  the  x  co-ordinate.  Since  the  solutions  obtained  for  each  case  will  in  general  be 
different,  the  arithmetic  mean  is  taken  as  an  approximate  solution.  There  are  of  course 
more  sophisticated  methods,  however  the  most  general  solution  appears  to  be  that 
derived  by  York[3].  The  York  algorithm  is  used  in  spg_fit,  and  is  derived  in  the  next 
section. 

A. 2  The  York  algorithm 

Let  xj  and  yj  be  the  actual  data  points  with  X}  and  Yj  denoting  the  fit  estimates,  i.e. 

Yi  =  a  +  bXi  (A.l) 


The  minimising  function  will  be  taken  as 


’  =  t  "  xi^  +  2 


(Vi  -  yi) 


: 


*i 


yi 


(A. 2) 


where  the  a  ,  and  a 

xi  Yi 


are  the  variances  for  each  point,  which  may  be  determined 


during  data  acquisition.  If  the  variances  are  not  known,  they  are  automatically  set  to 
unity  by  spg_fit. 


Following  the  usual  procedures,  minimising  equations  (A.l)  and  (A.2)  implies  the 
conditions 


5S  =  \  (Xi  '  Xi>5xi  +  \  (Yi  -  yi)5Yi  =  0 

i  ct  a 

xj  yi 


(A. 3) 


-5Yi  +  8a  +  bSXj  +  X[8b  =  0 


(A. 4) 
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York  then  multiplies  each  of  equations  (i=l,...)  (A.4)  by  its  own  Lagrangian 
undetermined  multiplier  Xj  and  takes  the  sum  of  the  ensuing  equations  plus 
equation  (A.3),  giving 


[bAj  +  _  (Xj  -  xj)]  SXj  +  [\i  +  2  (Yi  •  yi)]  5Yj 


*i 


yi 


+  8a /q  +  5b  )qXj  =  i 


(A.5) 


which,  by  equating  the  coefficients  to  zero,  implies  the  conditions 


Xi  —  "  Xjbcr  .  +  xj 
xi 

(A. 6a) 

Yi  =  Xi<7y.  +  yi 

(A. 6b) 

,  Ai  =  0 
i 

(A.6c) 

XXiXi  =0 

(A.6d) 

i 


Substituting  equations  (A.6a)  and  (A.6b)  into  equation  (A.l)  gives  the  expression  for  ?q 
as 


=  Wj  (a  +  bxj  -  yj)  (A .7) 

where 

wi-T/r  <A-8> 

b  V°yi 


Substituting  equation  (A .7)  into  equation  (A.6c)  gives  an  expression  for  the  intercept  as 


2,Wiyi-bX  WjXi 
i 

=  Y-bX 

where  the  "bar"  terms  are  defined  by 

XWiZi 


i 


l 


(A. 9a) 


(A.9b) 


(A.9c) 
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i.e.  the  line  will  go  through  the  "centre  of  gravity"  of  the  data. 

Introducing  the  quantities  U[  =  xj  -  X  and  Vj  =  yj  -  Y  ,  equation  (A.9b)  can  be  used  to 
simplify  equation  (A.7)  for  Xi  to 

Xi  =Wi(bUi-Vi)  (A. 10a) 

and  then  equation  (A.6a)  can  be  expressed  as 

Xi  =  -  Wi  (bUi  -Vi)  a*  +  xi  (A.IOb) 

X1 

Using  these  simplifications,  equations  (A.10)  are  substituted  into  equation  (A.6d)  and 
terms  in  powers  of  b  are  gathered.  The  result  is 

b3  +  3ab2  +  3[3b  +  y  =  0  (A. 11a) 


where 


2X  CTxi 

i 

“x,  Wi^i2 

i 


(A. lib) 


X  CTxi  Wi2Vi2-XUiWixi 

P-- - : - - -  (A. lie) 

3X  CTxi  Wi2ui2 


WiVixi 


Y  = 


X  <4  Wi2ui2 


(A. lid) 


Equation  (A.lla)  is  called  the  "least-squares  cubic"  by  York,  although  it  isn't  really  a 
cubic,  since  b  appears  in  the  definition  of  Wi-  The  form  of  equations  (A. 11c)  and  (A. lid) 
is  slightly  different  to  the  the  definitions  given  in  [3],  however  it  can  be  shown  that 
the  expressions  are  equivalent.  Equation  (A.lla)  can  be  solved  as  a  cubic  if  a  normal 
unweighted  least-squares  fit  is  performed  and  the  estimate  thus  obtained  for  b  used  in 
the  calculation  of  Wp  The  required  solution  to  equation  (A.lla)  is  given  by[2] 
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knew  =  -2(a2  -  p)1''2  cos(— ^)  -  a 


(A. 12a) 


where 


cos  <j>  = 


3  3  1 

“  ‘2a|3  +  2Y 

(a2  -  p)3/2 


(A. 12b) 


Equation  (A.12a)  is  the  third  root  [3]  of  the  least-squares  "cubic",  which  always 
appears  to  be  the  required  solution.  By  iterating  this  solution,  the  accuracy  can  be 
improved  to  any  desired  level.  The  convergence  properties  associated  with  the 
technique  of  treating  equation  (A.lla)  as  a  pseudo-cubic  and  iteratively  solving  for  b 
using  equation  (A.12a)  have  been  checked  for  a  number  of  sets  of  data  by  using  a  simple 
bisection  search  to  solve  for  the  appropriate  root  of  equation  (A.lla)  directly.  In  the 
experience  of  the  authors  the  pseudo-cubic  technique  has  been  reliable,  and  as  such  the 
technique  is  used  in  spg_fit.  Once  the  optimal  solution  is  achieved,  the  final  value  for 
b  is  then  used  one  last  time  to  calculate  the  corresponding  intercept,  a,  from 
equation  (A. 9a). 

Having  obtained  the  optimum  values  for  the  slope  and  intercept,  the  residuals  can  be 
calculated  from  equations  (A.6a),  (A.6b),  (A.7)  and  (A.IOa)  in  one  of  several  forms,  i.e. 

Xj  -  xi  =  bWj(bUj  -  Vi)  =  bWj(a  +  bxj  -  yj)  a2  (A. 13a) 

X1 


Yi  -  yi  =  Wi(bUi  -  Vi)  a*.  =  Wi(a  +  bxi  -  yi)  a2y. 


(A.13b) 


which  give  an  expression  for  the  minimisation  function  S  defined  by  equation  (A.2)  as 


=  XWi 


i"(bUi  -  Vi) 


2(u 2  2  2  ^ 
'b  a  +ct 
xi  yi 


-I 


Wiz(a  +  bxi  -  yj) 


2(u 2  2  2  ^ 

b  o  ,  n 
xj+  yj 

V  J 


(A. 14) 


The  program  spg_fit  displays  the  result  of  calculating  equation  (A.14)  divided  by  the 
degrees  of  freedom,  next  to  the  message  "reduced  chi  squared:”  which  appears  in  the 
"linear  fit"  panel  during  execution. 
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In  order  to  estimate  the  errors  in  the  fit  parameters,  York  uses  a  Taylor's  series 
expansion  for  the  straight  line  about  the  assumed  values  for  slope,  intercept  and 
estimated  positions  of  the  points[3].  The  resulting  expressions  are: 


Wi(bUi  -Vi)2 


2  1 
ab"n-2 


i 


Xwiui2 


(A. 15a) 


XWixi2 

i 

XWi 

i 


(A.15b) 


The  95%  confidence  error  bars  are  then  given  by  8a  =  2aa  and  8b  =  2 as  usual. 
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LEAST-SQUARES  MINIMISATION  OF  A  NON-LINEAR  FUNCTION 
B.l  Introduction 

In  this  section  an  outline  of  the  development  of  the  Levenberg-Marquardt  algorithm  for 
non-linear  least-squares  optimisation  and  the  method  of  determining  the  error 
estimates  for  the  fit  parameters  is  presented. 

B.2  The  x2  hypersurface 

Consider  an  arbitrary  function  y(x;a)  defined  at  the  points  xj  (i=l..N)  which  depends  in 
a  non-linear  manner  on  parameters  aj  (j=l..M).  The  parameters  a  are  to  be  chosen  so 
that  in  some  sense  they  model  the  set  of  experimental  data  points  yj  that  represent  the 
function  y(x;a),  each  of  which  has  a  measured  standard  deviation  aj.  In  order  to  model 
the  data,  it  is  desirable  to  minimise  the  usual  least-squares  norm 


[yi  -  y(xi;a)] 


2 


(B.l) 


as  a  function  of  the  parameters  a,  i.e 


dy2  a  Y  fyi  -  y(xi;a)]2 

5aj  3a;  2 

]  ]  i=l  CTj 


(B.2) 


The  quantity  j}  describes  an  M-dimensional  hypersurface  which  must  be  searched  to 
find  the  appropriate  minimum  value  of  x  •  Carrying  out  the  differentiation  yields 


N 

3X2  _V  [yj  -  y(xj;a)]  3y(xj;a) 

3a;  «  2  3ai 

1  i=l  a.  > 


(B.3) 


It  will  be  useful  later  to  have  the  matrix  of  second  derivatives  of  the  hypersurface  - 
the  Hessian  matrix.  This  is  given  by 
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The  term  [yi  -  y(xj;a)J  should  reflect  the  (hopefully)  small  random  measurement  error 
at  each  data  point,  and  as  such  should  not  contribute  significantly  to  the  summation  for 
good  models  and  reasonable  data.  The  effect  of  neglecting  the  second  derivative  term  is 
simply  to  alter  the  search  path  taken  along  the  hypersurface  to  the  x2  minimum  -  the 
value  of  the  minimum  itself  is  unaltered.  A  further  advantage  in  ignoring  the  last  term 
is  that  the  routine  is  made  more  robust  to  outlier  points. 


It  is  usual  to  define  from  equations  (B.3)  and  (B.4)  the  quantities 


Pj 


l^X2 
2  3ai 


(B.5) 


1 

a)k  =23aj3ak 


(B.6) 


The  symmetric  matrix  a  is  called  the  curvature  matrix  since  it  describes  the  curvature 
of  the  %2  hypersurface  in  parameter  space. 


Having  defined  the  function  to  be  minimised,  it  is  now  necessary  to  develop  a  search 
strategy  to  optimally  update  the  initial  guess  for  the  parameters  a. 

B.3  Finding  the  minimum 

By  definition,  the  gradient  Vx^,  or  that  direction  in  which  x2  increases  most  rapidly,  is 
a  vector  whose  components  are  equal  to  the  rate  at  which  x2  increases  in  that  direction. 
The  steepest  descent  search  algorithm  is  based  simply  upon  taking  small  steps  in  the 
opposite  direction  to  V%2  jn  parameter  space,  i.e. 

Y  9 

anew  =  Current  ^  ^X  'Current)  (B.7) 


where  anew  is  the  new  vector  of  parameters,  acurrent  is  the  vector  of  current  parameter 
Y 

values  and  ^  is  some  small  constant  which  determines  the  step  size.  Defining  the 
parameter  update  vector 

Sa  =  anew  *  Current  ( B .  8 ) 


then  equation  (B.7)  can  be  written  simply  as 


8a  =  yP 


(B.9) 


where  the  vector  (3  is  defined  by  equation  (4.5). 

It  is  well  known  however  that  the  steepest  descent  method  is  only  useful  when  the 
initial  parameter  vector  a  is  far  from  the  optimal  parameters  a0pt  that  minimise  x2- 
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When  a  =  a0pt,  the  steepest  descent  method  has  poor  convergence  properties. 
A  preferable  technique  is  to  assume  that  near  the  y2  minimum  the  hypersurface  is 
approximately  parabolic  with  respect  to  a  variation  in  the  parameters  a.  This  is 
equivalent  to  taking  the  first  three  terms  of  a  Taylor’s  series  expansion  about  some 

origin  a  near  the  minimum.  Any  other  point  nearby  in  parameter  space  can  then  be 
written  as 

A*>-  X2W  +  X  tp 

j  J  jk  J 

=  x^(a)-b*a  +  ^a»A*a  (B.10) 


where  b  =  | -  and  Ajk  =  ^^|-. 

Taking  the  gradient  of  equation  (B.10)  with  respect  to  the  parameters  aj  yields 
Vx2(a)  =  A»a  -  b  (B.ll) 


If  a  represents  the  optimal  set  of  parameters,  then  V%2  =  o  and 

A*aopt  =  b  (B.12) 

Subtracting  equation  (B.ll)  from  equation  (B.12)  and  re-arranging  gives  the  strategy  for 
determining  the  optimal  parameters  from  the  initial  vector  as 

aopt  -  a  -  A’l[Vx2(a)j  (B.13) 

This  technique  is  known  as  parabolic  extrapolation.  In  order  to  compare  this  result 
with  the  steepest  descent  algorithm,  equation  (B.13)  is  best  expressed  in  terms  of  the 
parameter  update  vector  8a,  the  vector  (3  and  the  matrix  a  as 


a.*8a  =  b 


or 


where  the  matrix  £  is  the  inverse  of  a. 


(B.14) 


(B.15) 


A  disadvantage  of  using  the  parabolic  extrapolation  technique  is  that  it  cannot  be 
relied  upon  to  approach  the  minimum  from  any  region  where  the  y}  hypersurface  is  not 
approximately  parabolic.  Specifically,  regions  where  the  curvature  is  negative  will 
lead  to  solutions  of  equation  (B.15)  that  are  unreliable. 
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The  choice  of  using  either  the  steepest  descent  or  the  parabolic  extrapolation  technique 
is  clearly  dependent  upon  the  accuracy  with  which  the  initial  conditions  can  be 
specified,  and  the  quality  of  the  experimental  data.  If  the  initial  values  of  the 
parameters  are  not  well  known,  then  the  steepest  descent  algorithm  is  the  more  robust 
method  of  finding  the  minimum.  Once  the  minimum  is  approached  however,  use  of 
the  parabolic  extrapolation  technique  is  preferred. 


B.4  The  Levenberg-Marquardt  algorithm 

The  Levenberg-Marquardt  (or  gradient  expansion)  algorithm  combines  the  best  features 
of  both  steepest  descent  and  parabolic  extrapolation  techniques.  Recalling  the 
parameter  update  formulae  for  each  method 


(steepest  descent) 


oc»8a  =  p 


(parabolic  extrapolation) 


then  a  new  matrix  a  is  constructed  by  defining 


5jpajj(l+W  (B.16) 

ajk^ajk  ;  j*k  (B.17) 

where  A.  is  a  scalar  quantity  known  as  the  Marquardt  parameter.  The  new  parameter 
update  formula  now  takes  the  form 


a«Sa  =  p  (B.18) 

When  X  is  very  small  or  zero,  equation  (B.18)  is  just  the  parabolic  extrapolation 
technique.  However,  when  X  is  large  compared  to  the  oijk,  the  matrix  a  is  forced  to  be 
diagonally  dominant,  and  equation  (B.18)  becomes 

ajjX5a  =  p  (B.19) 

1 

which  is  just  the  steepest  descent  algorithm  with  ctjjX.  <=>  — .  Thus,  by  adjusting  the  size 

of  X,  the  technique  can  be  forced  into  either  a  steepest  descent  or  a  parabolic  search  of 
the  hypersurface.  The  formal  solution  for  the  parameter  update  vector  is  then 

8a  =  P»2  (B.20) 


where  £  is  the  inverse  of  the  pseudo-curvature  matrix  a. 
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The  Levenberg-Marquardt  algorithm  proceeds  in  the  following  manner: 

(1)  calculate  %^(a) 

(2)  start  with  a  small  value  of  X,  typically  0.001 

(3)  solve  equation  (B.20)  for  8a  and  evaluate  x^Ca+Sa) 

(4)  if  %2(a+sa)  >  ^2(a)  increase  X  by  a  substantial  factor  (typically  10)  and  go  to 
step  (3) 

(5)  if  x^(a+8a)  <  x^(a)  decrease  X  by  a  substantial  factor  (typically  10),  update 
the  trial  solution  a+Sa  — »  a,  and  go  to  step  (3) 


The  routine  can  be  terminated  when  all  parameters  have  sufficiently  converged,  when 
does  not  significantly  decrease  or  after  a  maximum  number  of  iterations.  A  further 
check  for  good  convergence  is  that  the  Marquardt  parameter  be  small  compared  to 
the  ajj. 

B.5  Error  estimation  in  non-linear  least-squares 

Since  the  result  of  the  least-squares  minimisation  of  a  non-linear  function  requires  some 
form  of  search  routine  rather  than  an  exact  analytic  solution,  there  is  no  analytical 
form  for  the  estimated  error  in  each  of  the  fit  parameters  aj.  The  best  that  can  be 
achieved  is  to  develop  an  algorithm  which  is  reasonable  and  reduces  to  the  correct 
expression  for  a  function  that  is  linear  in  a. 


In  order  to  achieve  this,  it  is  appropriate  to  briefly  review  the  derivation  of  the  exact 
analytic  expression  for  the  estimated  errors  as  they  occur  in  the  case  of  linear  least- 
squares.  In  this  case  the  optimal  parameters  a  are  determined  from  the  linear 
equivalent  of  equation  (B.15),  which  is 


aj  =  £jk  Pk 

k 

where  the  model  is 

y(*i)  =  ^  aj  bj(xi) 


with 


N 


Pk^ 


ld%2  Vyi 

2  3a]< 


bk(xi) 


i=l  ct. 


and 


o  i  •  u  1  92X2 

8  =  a"1  with  ajk  =  -^-Trr 


2  3aj3ak 


(B.21) 


(B.22) 


(B.23a) 
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N 


i=l  a. 


'  bj(xi)  bk(xi) 


(B.23b) 


Now,  the  standard  deviation  a  representing  the  uncertainty  in  the  determination  of 


any  parameter  aj  is  the  square  root  of  the  sum  of  the  squares  of  the  product  of  the 
standard  deviation  of  each  data  point  cq  multiplied  by  the  effect  that  data  point  has 
on  the  determination  of  parameter  aj[l],  i.e. 


(B.24) 


From  equation  (B.21) 


(B.25) 


since,  from  equation  (B.23b),  ajk  is  independent  of  yi,  and  so  8  =  a'1  is  independent  of  yi 
also.  From  equation  (B.23a) 


aPk 

3yi 


rbk(xi) 


(B.26) 


and  so 


M 

8a;  V  1 

1  a2  bk£jk 

k=l  CTj 


(B.27) 


Substituting  equation  (B.27)  into  equation  (B.24)  yields 
MM  f  N 

naj  Z 

1  k=l  r=l  ^i=l 

M  M 

£jk  £jr  «kr  (by  equation  (B.23b)) 

k=l  r=l 


A 


bk(xj)  br(xi) 


ai 


J 
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M 

£jk  Sjk 
k=l 


i.e. 


(B.28) 


Thus,  in  the  linear  case,  the  diagonal  elements  of  the  inverse  curvature  matrix  are  the 
variances  of  the  fitted  parameters  a.  The  off  diagonal  elements  £jk  (k*j)  describe  the 
covariances  between  the  parameters  aj  and  afc.  The  matrix  S  is  therefore  known  as  the 
covariance  matrix  or  error  matrix.. 


Comparison  of  the  linear  and  non-linear  solutions  (equations  (B.21)  and  (B.15))  shows 
that  the  fundamental  mathematical  difference  between  the  two  problems  is  that  in  the 
non-linear  case  the  parameter  increments  are  calculated,  rather  than  the  parameters 
themselves.  The  error  in  the  linear  case  is  associated  solely  with  the  curvature  matrix 
a,  and  from  equations  (B.23b)  and  (B.6)  it  is  apparent  that  the  definition  of  a  is 
identical  for  both  problems. 

The  approximation  made  in  the  non-linear  case  is  that  the  curvature  matrix  calculated 
at  the  x2  minimum  still  reflects  the  variances  of  each  of  the  fitted  parameters,  and  so 
the  final  step  in  the  Levenberg-Marquardt  routine,  once  satisfactory  convergence  is 
achieved,  is  to  set  the  Marquardt  parameter  X  equal  to  zero  and  to  re-calculate  a  and 
then  8.  The  diagonal  elements  £jj  are  taken  to  be  the  estimates  of  the  variance  for  each 
parameter  aj. 

2 

Equation  (B.28)  assumes  that  the  variance  cr  of  each  data  point  is  known.  If  they  are 
2 

not  known  (i.e.  cr  =1),  they  must  be  estimated  from  the  experimental  data  set  using 

N 

CTi=s2  =  N^lX  [yi  -  y(xi;a)]2 
i=l 


1  2 

"  N-MZ(1) 


(B.29) 
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where  %q)  denotes  the  sum  of  the  squares  of  the  residuals  calculated  with  cq  =  1.  The 
variance  of  the  parameter  aj  is  then  calculated  from 

2 

2  9  *(1) 

aa.  =  s2ejj(ai=1)=— Ejjfa^D  (B.30) 

The  95%  confidence  interval  for  each  parameter  is  then  just 


5aj  =  2a. 
j  aJ 


(B.31) 


irrespective  of  whether  equation  (B.28)  or  equation  (B.30)  is  used  to  calculate  a. 


a)' 
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/* 


*/ 


APPENDIX  C 

THE  README  AND  spg_fit  HEADER  FILES 

SPG-Fit 

Copyright  DSTO,  1990 

Authors:  Shaun  Troedson  Tony  Lindsay 
SP  Group 

Communications  Division 
DSTO  Salisbury 


This  code  may  be  copied,  modified  etc  provided  the  above  copyright  and 
authors'  names  remain  above  the  major  routines.  Any  bugs  or  suggestions 
may  be  sent  to  the  above  address  or  try  email  on  sct@dstos3 . dsto . oz 

Note  that  the  source  files  should  be  placed  in  a  general  access  read 
write  directory  since  there  will  be  frequent  recompiling  by  the  user. 


Setting  up  the  environment: 

To  set  up  the  environment  do  the  following: 

1.  Create  a  source  file  containing  the  model  function 
or  just  copy  from  the  functions  provided  in  this 
directory.  This  is  not  essential. 

2 .  Optionally  include  in  the  same  file  the  1st  derivatives  of 
this  model  w.r.t  the  parameters.  This  will  improve  speed 
and  reduce  round  off  errors. 

3.  Create  a  text  file  (fn_info)  with  the  parameter  names. 

The  provided  functions  have  info  files  in  this  directory 
with  the  prefix  (info_) . 

4.  Recompile  using  a  simple  make  command. 

Step  1,2: 

The  user  model  may  be  written  in  Pascal  or  C.  The  filenames 

are  "cfn.c"  or  "pfn.p"  for  the  C  or  Pascal  routine  respectively. 

The  main  procedure  should  be  called  fn (N, X, P, Y, deriv, dY)  where 

N  is  the  number  of  data  points; 

X  is  an  array  of  size  NPTS  (see  spg_fit.h)  of  real  (double) 
that  contains  the  independent  variable  values. 

If  the  user  fn  is  not  in  closed  form,  ie  it  is 
not  possible  to  simply  plug  these  values  in  to 
get  the  function  value  at  that  point,  care  should 
be  taken  that  the  values  for  the  dependent  variable 
match  the  corresponding  values  for  X; 

P  is  an  array  of  size  NPARAMS  (see  spg_fit.h)  of  real  (double) 
that  contains  the  values  of  the  parameters; 

Y  is  an  array  of  size  NPTS  of  real  (double)  that  contains 
the  theoretical  curve  dependent  variable  values  and 
is  set  by  the  user  program; 

deriv  is  a  flag  which  is  only  optionally  used  by  the  user 
program  to  determine  whether  to  perform  analytic  derivative 
calculations.  spg_fit  calls  the  user  routine  to 
obtain  values  to  plot  and  does  not  require  the 
derivative  information,  hence  in  order  to  save 
time  the  user  should  use  this  parameter  to 
conditionally  perform  the  derivative  calculation 
when  this  flag  is  1. 

dY  is  a  matrix  of  size  [NPARAMS, NPTS]  of  real  (double) 
returned  by  the  user  program.  Each  row  contains  the 
derivatives  of  the  of  the  function  w.r.t  each  parameter 
over  the  values  of  the  array  X.  This  is  optional  and 
if  spg_fit  does  not  detect  that  the  analytic  derivatives 
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exist  it  will  perform  a  numeric  calculation.  It  should  be 
noted  however  that  this  entails  taking  the  difference 
of  two  similar  numbers  and  hence  may  suffer  due  to 
round  off  error. 


Step  3 : 

The  file  "fn_info"  is  a  text  file  that  contains  the  parameter 
names-  line  by  line.  Note  that  these  should  be  in  the  same  order 
as  implied  in  the  array  P  above.  If  the  supplied  routines  are 
used  then  their  corresponding  info  file  can  be  used. 

Step  4:  Compiling. 

The  fitting  routine  is  now  ready  for  compiling.  The  file  Makefile 
contains  two  variables  which  should  be  defined.  RESIDE  is  the 
directory  in  which  the  source  codes  and  object  files  for  spg_fit  live. 
DEST  is  the  directory  in  which  the  user  wishes  to  place  the 
final  executable.  The  file  "spg_fit.h"  contains  other  definitions 
which  may  be  specified.  Here  the  maximum  number  of  points  (NPTS) 
and  maximum  number  of  parameters  (NP ARAMS)  are  defined.  Also 
defined  are  STEPSIZE  (the  step  taken  in  calculating  the  numerical 
derivatives) ,  CRITERION  (the  convergence  criterion) ,  and  the 
local  fonts  library. 

Finally  a  simple  make  command  is  used  to  compile  and  link. 

If  you  are  in  the  source  directory  for  spg_fit  then  type 

make  fitp  — or —  make  fitc 

if  the  user  routine  is  in  Pascal  or  C  respectively.  You  may 
also  use  the  -f  option  of  make  if  you  are  not  in  the  source 
directory.  The  routine  may  now  be  executed  by  typing 

fit 


The  Data  file: 

The  format  of  the  text  file  consists  of  2,3  or  4  columns.  The  first 
column  is  taken  as  the  independent  variable  values  and  the  second 
column  the  dependent  variable  data.  The  3rd  and  4th  columns,  which 
are  optional,  contain  the  variances  of  the  individual  data  points. 
These  columns  are  used  when  performing  a  weighted  fit, though  the  4th 
column  is  only  used  when  performing  a  X  &  Y  weighted  linear  fit. 

The  3rd  column  is  taken  to  be  the  variances  in  the  independent 
column,  or,  if  the  4th  column  does  not  exist,  the  variances  in  the 
dependent  column.  The  columns  may  be  separated  by  spaces  or  tabs. 
Note  that  the  data  should  be  presorted  otherwise  the  results  are 
unpredictable . 

Some  notes  on  operation: 

The  routine  allows  parameters  to  be  changed  on  the  screen. 

However  it  only  updates  the  chi  squared  value  if  the  number 
of  samples  is  equal  to  the  number  of  data  points  (this  is  the 
default) .  This  is  so  that  two  calls  to  the  function  are  not 
required  which  may  be  annoying  for  slow  functions. 

The  routine  uses  a  Levenberg-Marquardt  routine  for  the  non¬ 
linear  fitting,  hence  the  field  for  lambda.  The  linear  fitting 
uses  a  York  routine  to  do  a  X  &  Y  weighted  linear  fit. 

References  may  be  found  in  the  documentation. 
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Mar  8  10:03  1991  spg_fit.h  Page  1 


/**  NPTS :  maximum  number  of  allowed  data  and  plot  points 

hopefully  this  will  be  malloc'd  in  a  later  version*/ 
♦define  NPTS  250 

/**  NP ARAMS :  maximum  number  of  allowed  paramters*/ 

♦define  NP ARAMS  10 

/**  STEPSIZE:  this  is  the  divisor  of  a  parameter  in  calculating 
the  derivatives  of  chi  sq  wrt  this  parameter.  This 
constant  may  be  optimized  for  the  application*/ 

♦define  STEPSIZE  le6 

/ * *CRITERION :  the  convergence  criterion,  fitting  will  stop  if  the 
fractional  change  of  all  parameters  between  one  iteration 
and  the  next  is  less  than  this  number*/ 

♦define  CRITERION  le-5 

/**  NSD:  number  of  standard  deviations  to  use  for  errors*/ 

♦define  NSD  2 


/ **MAX_ITERATIONS :  if  the  LM  routine  gets  into  a  runaway  process 
it  will  stop  at  this  value*/ 

♦define  MAX_ITERATIONS  50 

/**FONTS:  this  defines  your  font  directory.  BOLD_FONT  will  be  the 
font  used  for  field  labels.  PLAIN_FONT  will  be  the 
font  used  for  field  values.*/ 

♦define  BOLD_FONT  " /usr/ lib/f onts/f ixedwidthf onts/cour . b . 12 " 

♦define  PLAIN  FONT  "/usr/lib/fonts/f ixedwidthfonts/cour . r . 12" 
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APPENDIX  D 


A  LISTING  OF  SUPPLIED  FUNCTIONS  INCLUDING  PASCAL 


AND  C  ROUTINES 

Mar  8  10:34  1991  f it_gauss_area .p  Page  1 


/* 

gaussian  distribution  function 
total  area  normalized  to  P[4] 

FWHM  =  P  [2] 
centred  at  P[l] 
baseline  of  P[3]; 

*/ 

const  NPTS=250;  /*take  care  that  these  values  are  the  same  as*/ 

NPARAMS=10;  /*in  the  spg_fit.h  file*/ 

f ourln2=2 .7725887; 
sqrtpi=l . 772453851; 

type 

parray=array [1 . .NP ARAMS]  of  real;  /*parameter  type  array*/ 

darray=array [1 . .NPTS]  of  real;  /*data  points  type  array*/ 

derivarray=array [1 . .NP ARAMS,  1 . .NPTS]  of  real;  /*derivative  type  array*/ 


procedure  f n (N : integer ; 

var  x:darray; 
var  P : par ray; 
var  F:darray; 
deriv: integer; 
var  dF:derivarray)  ; 


/‘number  of  data  points*/ 

/♦independent  values*/ 

/♦parameter  values*/ 

/♦dependent  values  (returned) */ 

/♦flag  for  derivatives  to  be  performed*/ 
/♦derivative  values  (returned) */ 


var 

i : integer; 
temp : da r ray; 
u : real; 

FWHM: real; 
base : real; 
area : real; 

begin 


u :  =P  [  1 ] ; 
FWHM: =P [2] ; 
base :=P  [3]  ; 
area:=P  [4]  ; 


/♦mean*/ 

/♦full  width  half  maximum*/ 
/♦offset  from  zero*/ 

/♦area  under  curve*/ 


for  i:=l  to  N  do  /*evaluate  function*/ 

begin 

temp  [i] :=sqrt (fourln2) *exp (-fourln2*sqr ( (x [i] -u) /FWHM) ) / sqrtpi ; 
F ( i ] :=base+area* (temp [i] /FWHM) ; 

end; 


if  deriv=l  then  /*optional*/ 

for  i:=l  to  N  do  /*evaluate  derivative  w.r.t  parameters*/ 

begin 

dF [ 1] [i] : =2*f ourln2* (x [i] -u) *area*temp [ i] / (FWHM*sqr (FWHM) ) ; 
dF [2] [i] := (2*f ourln2*sqr ( (x [i] -u) /FWHM) -1) *area*temp [i] /sqr 
dF [3] [i] :=1 . 0; 
dF [4] [i] :=temp[i] /FWHM; 

end; 


end; 
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10:29  1991  fit  lorentz_area . c  Page  1 


/* 

lorentzian  distribution  function 

total  area  under  curve  normalized  to  P [ 3 ] 

P [0]=FWHM; 

centred  at  P[l] 

baseline  P[2]; 

*/ 


finclude  <math.h> 
tinclude  "spg_fit.h" 

fn (n,  X, P, f , deriv, df ) 
int  n; 

double  X [ ] ; 
double  P [ ] ; 
double  f [ ] ; 
int  deriv; 

double  df [] [NPTS] ; 

{ 

int  i  ; 

double  temp [NPTS]; 
double  Pi=3 . 1415926538; 

for  ( i=0 ; i<n; i++)  /*evaluate  function*/ 

temp [i]=4.0*(X[i]-P[l])*(X[i]-P[l])+P[0]*P[0]; 
f [i]=P [2] +2 ,0‘P [3] *P [0] / (Pi*temp [i] ) ; 

} 

if  (deriv)  /*optional*/ 

for  (i=0;i<n;i++)  /*evalaute  derivatives  w.r.t  params*/ 

df  [0]  [ i]  =  ( 1 .0-2 .  0*P [0] *P [0] /temp [i] ) *2.0*P [3] / (Pi*temp[i]  ) 
df  [1]  [i]=16*P [3] *P [0] * (X[i]-P [1] ) / (Pi*temp [ i] *temp[i] ) ; 
df [2] [i] =1 . 0 ; 

df [3] [i]=2.0*P(0]/ (Pi* temp [i] ) ; 


/*spg_fit  header  file*/ 


/‘number  of  data  points*/ 

/‘independent  values*/ 

/‘parameter  values*/ 

/‘dependent  values  (returned) */ 

/‘flag  for  derivatives  to  be  performed*/ 
/‘derivative  values  (returned)*/ 
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Mar  8  10:31  1991  fit_exp.c  Page  1 


/* 

exponential  function:  y[i]  =  offset  +  Y [0 ] *exp (exponent *x [i]  )  ; 

*/ 


♦include  <math.h> 
♦include  "spg_fit.h" 

fn (n, X,  P, f , deriv, df ) 
int  n  ; 

double  X [ ] ; 
double  P  [  ]  ; 
double  f  [  ]  ; 
int  deriv; 

double  df[][NPTS]; 


/‘number  of  data  points*/ 

/* independent  values*/ 

/♦parameter  values*/ 

/♦dependent  values  (returned) */ 

/♦flag  for  derivatives  to  be  performed*/ 
/♦derivative  values  (returned)*/ 


{ 


int  i  ; 


I 


for  (i=0 ; i<n; i++)  /‘evaluate  function*/ 

f  [i]=P[2]+P[0]*exp(P[l]*X[i])  ; 


if  (deriv)  /‘optional*/ 

for  (i=0; i<n; i++)  /‘evaluate  derivs  w.r.t  params*/ 

{ 

df  [0]  [i]=exp (P [1] *X[i] ) ; 

df  [1]  [i]=P [0] *X[i] *exp (P [1] *X[i] ) ; 

df [2]  [i] =1 . 0 ; 

} 
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Mar  8  10:33  1991  f it_gauss_ht .p  Page  1 


/* 

gaussian  distribution  function 
normalized  to  height  P[4] 

P [2 ]  is  the  FWHM 
centred  at  P [ 1 ] 
baseline  P  [3] 

*/ 

const  NPTS=200;  /*take  care  that  these  values  are  the  same  as*/ 

NPARAMS=10;  /*in  the  spg_fit.h  file*/ 

fourln2=2 .7725887; 


type 


parray=array [1 . .NP ARAMS]  of  r 
darray=array [1 . .NPTS]  of  real 
derivarray=array [1 . .NP ARAMS, 1 

procedure  f n (N : integer ; 

var  x:darray; 
var  P:parray; 
var  F:darray; 
deriv : integer; 
var  dF :derivarray)  ; 


al;  /*parameter  type  array*/ 

/*data  points  type  array*/ 

.NPTS]  of  real;  /*derivative  type  array*/ 

/‘number  of  data  points*/ 

/‘independent  values*/ 

/‘parameter  values*/ 

/‘dependent  values  (returned) */ 

/‘flag  for  derivatives  to  be  performed*/ 
/‘derivative  values  (returned) */ 


var 

i : integer; 
temp:darray; 
u: real; 

FWHM: real; 
base : real; 
height : real; 

begin 

u : =P  [  1 ]  ; 

FWHM : =P  [2] ; 
base :=P  [3] ; 
height : =P [ 4] ; 

for  i:=l  to  N  do  /‘evalaute  function*/ 

begin 

temp [i] :=exp(-fourln2*sqr ( (x[i]-u) /FWHM) ) ; 

F[i] :=base+height*temp [i] ; 

end; 

if  deriv=l  then  /‘optional*/ 

for  i:=l  to  N  do  /‘evaluate  derivs  w.r.t  params*/ 

begin 

dF [1] [i] :=2*fourln2* (x[i] -u) ‘height *temp [i] /sqr (FWHM) ; 
dF [2]  [i] : =2 *f ourln2*sqr (x [i] -u) ‘height* temp  f i] / (FWHM* sqr (FWHI 
dF [3] [i] :=1 . 0; 
dF  [4]  [i]  :=temp [ i] ; 

end; 

end; 
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Mar  8  10:35  1991  f it_lorentz_ht . c  Page  1 


/* 

lorentzian  distribution  function 
normalized  to  P [ 3 ]  at  line  centre 
FWHM  =  P [ 0 ] 

centred  at  P [ 1 ]  with  baseline  of  height  P [ 2 ] 

*/ 


♦include 

<math . h> 

♦include 

"spg  fit . h" 

fn (n,  X,  P 

,  f, deriv,  df ) 

int 

n; 

/‘number  of  data  points*/ 

double 

X[]  ; 

/‘independent  values*/ 

double 

p  []  ; 

/‘parameter  values*/ 

double 

f  []; 

/‘dependent  values  (returned) */ 

int 

deriv; 

/‘flag  for  derivatives  to  be  performed*/ 

double 

df []  [NPTS ]  ; 

/‘derivative  values  (returned) */ 

{ 

int  i  ; 

double  temp [NPTS]; 
double  Pi=3 . 1415926538 

r 

for  ( i=0 ; i<n; i++) 

/‘evaluate  function*/ 

{ 

temp [i] =4 . 0* (X[i]-P [1] ) * (X[i]-P [1] ) +P [0] *P [0]  ; 
f  [i]=P [2]+P [3] *P [0] *P [0] /tempfi] ; 

} 

if  (deriv)  /*optional*/ 

for  (i=0; i<n; i++)  /‘evaluate  derivs  w.r.t  params*/ 

{ 

df [0] [ i ] = ( 1 . 0-P [0] *P [0] /temp[i] ) * (2*P [3] *P [ 0 ] /temp [ i] ) ; 
df [1] [ i] =8* (X[i]-P [1] ) *P [3] *P [0] *P [0] / (temp [ i] *temp [ i] ) 
df [2] [ i] =1 . 0 ; 

df [3] [i]=P[0]*P[0] /temp [i] ; 
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Mar  8  10:37  1991  f it_power_fn . c  Page  1 


/* 

power  function:  y[i]  =  intercept+scale*x [i] ''exponent 

*/ 

#include  <math.h> 

♦include  ”spg_fit.h" 

f n (n, X, P , f , deriv, df ) 
int  n; 

double  X [ ] ; 
double  P [ ] ; 
double  f  [  ] ; 
int  deriv; 

double  df [] [NPTS] ; 

{ 

int  i  ; 

double  powx[NPTS]; 
double  lnx[NPTS]; 

for  (i=0;i<n;i++)  /*evaluate  function*/ 

{ 

lnx [i] =log (X [i] ) ; 

powx [ i ] =exp ( P [ 0 ] * lnx [ i ] ) ; 

f [ i ] =P [ 2 ] +P [ 1 ] *powx [ i ] ; 

} 

if  (deriv)  /‘optional’*/ 

t 

for  (i=0; i<n; i++)  /‘evaluate  derivs  w.r.t  params*/ 

{ 

df  [0]  [i]=P [1] *lnx[i] *powx[i] ; 
df  [1]  [i]=powx[i] ; 
df  [2]  li]— 1.0; 

} 

} 

) 


/‘number  of  data  points*/ 

/‘independent  values*/ 

/‘parameter  values*/ 

/‘dependent  values  (returned)*/ 

/‘flag  for  derivatives  to  be  performed*/ 
/‘derivative  values  (returned)*/ 
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/* 

general  sinusoid:  y[i]  =  intercept  +  scale  *  sin(  k*x[i]  +  phase  ) 

*/ 

tinclude  <math.h> 
finclude  "spg_fit.h" 

fn (n,  X, P, f, deriv, df ) 
int  n; 

double  X [ ] ; 
double  P  [  ]  ; 
double  f  [  ]  ; 
int  deriv; 

double  df [] [NPTS] ; 

{ 

int  i  ; 

for  (i=0 ; i<n; i++)  /‘evaluate  function*/ 

f [i]=P [2]+P [0] *sin(P [1] *X[i]+P [3] ) ; 

if  (deriv)  /*optional*/ 

{ 

for  (i=0 ; i<n; i++)  /‘evaluate  derivs  w.r.t  params*/ 

{ 

df  [0]  [ i] =sin (P [ 1] *X[i]+P[3] ) ; 

df [1] [i]=P[0]*x(i]*cos(P(l]*Xti]+P[3]) ; 

df [2] [i] =1 . 0 ; 

df  [3]  [i]=P(0]*cos(P[l]*X[i]+P[3])  ; 

) 

} 

} 


/*spg_fit  header  file*/ 


/★number  of  data  points*/ 

/★independent  values*/ 

/★parameter  values*/ 

/★dependent  values  (returned) */ 

/★flag  for  derivatives  to  be  performed*/ 
/★derivative  values  (returned) */ 
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